/**
   parts.cu - Created by Timothy Morey on 8/15/2011.

 */

#include <cmath>
#include <cstdlib>
#include <ctime>
#include <stdio.h>
#include <string>

using namespace std;

#include "parts.cuh"

__global__ 
void UpdatePositionAndVelocity(DoF* X, DoF* Y, 
                               DoF* VX, DoF* VY,
                               DoF* AX, DoF* AY, 
                               unsigned int n, float dt,
                               float minX, float maxX,
                               float minY, float maxY,
			       unsigned int minPart, unsigned int maxPart)
{
  int i = blockIdx.x * blockDim.x + threadIdx.x;

  if(i + minPart <= maxPart)
  {
    DoF x = X[minPart + i];
    DoF y = Y[minPart + i];
    DoF vx = VX[i];
    DoF vy = VY[i];
    DoF ax = AX[i];
    DoF ay = AY[i];

    DoF dx = vx * dt + 0.5f * ax * dt * dt;
    DoF dy = vy * dt + 0.5f * ay * dt * dt;
    x += dx;
    y += dy;
    EnsureInBounds(x, y, minX, maxX, minY, maxY);

    vx += 0.5f * ax * dt;
    vy += 0.5f * ay * dt;

    X[minPart + i] = x;
    Y[minPart + i] = y;
    VX[i] = vx;
    VY[i] = vy;
  }
}

__global__
void UpdateAcceleration(DoF* X, DoF* Y, 
                        DoF* AX, DoF* AY, 
                        unsigned int n,
                        float minX, float maxX,
                        float minY, float maxY,
			unsigned int minPart, unsigned int maxPart)
{
  int i = minPart + blockIdx.x;
  int j = threadIdx.x;
  int otherPart = j;

  __shared__ DoF FX[THREADS_PER_BLOCK];
  __shared__ DoF FY[THREADS_PER_BLOCK];

  FX[j] = 0.0f;
  FY[j] = 0.0f;

  if(i <= maxPart && j < n)
  {
    DoF xi = X[i];
    DoF yi = Y[i];

    while(otherPart < n)
    {
      if(i != otherPart)
      {
  	DoF dx = XDistanceTo(xi, X[otherPart], minX, maxX);
  	DoF dy = YDistanceTo(yi, Y[otherPart], minY, maxY);
	
  	DoF r2 = (dx * dx) + (dy * dy);
  	DoF rm2 = 1.0f / r2;
  	DoF rm6 = rm2 * rm2 * rm2;
  	DoF foverr = 24.0f * rm6 * (2.0f * rm6 - 1.0f) * rm2;

  	FX[j] += foverr * dx;
  	FY[j] += foverr * dy;
      }

      otherPart += blockDim.x;
    }
  }

  __syncthreads();
    
  int redux = THREADS_PER_BLOCK / 2;
  while(redux != 0)
  {
    if(j < redux && j + redux < n)
    {
      FX[j] += FX[j + redux];
      FY[j] += FY[j + redux];
    }
    
    __syncthreads();
    redux /= 2;
  }
  
  if(i < n)
  {
    AX[i-minPart] = FX[0];
    AY[i-minPart] = FY[0];
  }
}

__global__ 
void UpdateVelocity(DoF* VX, DoF* VY, 
                    DoF* AX, DoF* AY,
                    float dt,
                    unsigned int minPart, unsigned int maxPart)
{
  int i = blockIdx.x * blockDim.x + threadIdx.x;

  if(minPart + i <= maxPart)
  {
    VX[i] += 0.5f * AX[i] * dt;
    VY[i] += 0.5f * AY[i] * dt;
  }
}

__global__ 
void CalculateMagnitude(DoF* VX, DoF* VY, DoF* VM,
                        unsigned int minPart, unsigned int maxPart)
{
  int i = blockIdx.x * blockDim.x + threadIdx.x;

  if(minPart + i <= maxPart)
  {
    DoF vx = VX[i];
    DoF vy = VY[i];
    VM[i] = sqrt(vx * vx + vy * vy);
  }
}

__global__
void ClearFrame(unsigned char* frame, unsigned int width, unsigned int height,
                unsigned char r, unsigned char g, unsigned char b, unsigned char a)
{
  int row = blockIdx.x;
  int col = threadIdx.x;
  while(col < width)
  {
    int i = row * width + col;
    frame[i*4+0] = r;
    frame[i*4+1] = g;
    frame[i*4+2] = b;
    frame[i*4+3] = a;
    col += blockDim.x;
  }
}

__global__
void RenderFrame(unsigned char* frame, unsigned int width, unsigned int height,
                 DoF* X, DoF* Y, unsigned int n,
                 float minx, float maxx, float miny, float maxy)
{
  int i = blockIdx.x;
  DoF x = X[i];
  DoF y = Y[i];
  int col = (int)floor((x - minx) / (maxx - minx) * width);
  int row = (int)floor((y - miny) / (maxy - miny) * height);
  int fi = row * width + col;
  frame[fi*4] = 255;
  frame[(fi+1)*4] = 255;
  frame[(fi-1)*4] = 255;
  frame[(fi-width)*4] = 255;
  frame[(fi+width)*4] = 255;
}

__host__ 
__device__ 
DoF XDistanceTo(DoF x1, DoF x2, float minX, float maxX)
{
  DoF retval = x1 - x2;

  DoF width = maxX - minX;
  if(retval > (0.5f * width))
    retval -= width;
  else if(retval < (-0.5f * width))
    retval += width;

  return retval;
}

__host__ 
__device__ 
DoF YDistanceTo(DoF y1, DoF y2, float minY, float maxY)
{
  DoF retval = y1 - y2;

  DoF height = maxY - minY;
  if(retval > (0.5f * height))
    retval -= height;
  else if(retval < (-0.5f * height))
    retval += height;

  return retval;
}

__host__ 
__device__ 
DoF DistanceTo(DoF x1, DoF y1, DoF x2, DoF y2,
               float minX, float maxX, float minY, float maxY)
{
  DoF dx = XDistanceTo(x1, x2, minX, maxX);
  DoF dy = YDistanceTo(y1, y2, minY, maxY);
  return sqrt(dx * dx + dy * dy);
}

__host__
__device__
void EnsureInBounds(DoF& x, DoF& y, float minX, float maxX, float minY, float maxY)
{
  DoF width = maxX - minX;
  DoF height = maxY - minY;

  if(x > maxX)
    x -= ceil((x - maxX) / width) * width;
  
  if(x < minX)
    x += ceil((minX - x) / width) * width;
  
  if(y > maxY)
    y -= ceil((y - maxY) / height) * height;
  
  if(y < minY)
    y += ceil((minY - y) / height) * height;
}
